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Abstract. We revisit a totaUy asymmetric simple exclusion process (TASEP) with 
open boundaries and a global constraint on the total number of particles [Adams, et. 
al. 2008 J. Stat. Mech. P06009]. In this model, the entry rate of particles into the 
lattice depends on the number available in the reservoir. Thus, the total occupation 
on the lattice feeds back into its filling process. Although a simple domain wall theory 
provided reasonably good predictions for Monte Carlo simulation results for certain 
quantities, it did not account for the fluctuations of this feedback. We generalize 
the previous study and find dramatically improved predictions for, e.g., the density 
profile on the lattice and provide a better understanding of the phenomenon of " shock 
localization." 

E-mail: lacooklSvt . edu , rkpziaOvt . edu 

Keywords: non-equilibrium statistical physics, totally asymmetric exclusion process, 
biological transport 

1. Introduction 

The Totally Asymmetric Simple Exclusion Process (TASEP) is a simple, yet rich 
model in the poorly understood and vast realm of non-equilibrium statistical mechanics. 
Particles are placed in a hypercubic lattice (for example) and hop randomly to a nearest 
neighbor empty site, except along one of the axes where the hopping is uni- directional. 
Since its dynamics violate detailed balance, its stationary states are non-trivial, with 
typically no equivalence to any states in equilibrium statistical mechanics. For a system 
with periodic boundary conditions, the stationary state is trivial, with all configurations 
having equal probability (i.e., a flat distribution) . However, the dynamic properties 
are far from trivial, being quite distinct from those for a completely symmetric exclusion 
processes (i.e., simple diffusion) [2] . For systems with open boundaries, in which 
particles hop in and out of the system with various rates, even the stationary states 
are quite complex. Indeed, an open TASEP in just one dimension settles into three 
different phases, depending on the entry and exit rates[3J. Despite the simplicity 
of its microscopic rules of evolution, the solution to this process was not known 
analytically until recently [H El E]. Needless to say, its dynamic properties are even 
more complex^ [HI [9l [10]. Examples of comprehensive reviews on this "simple" process 
include pTl [T2] . Meanwhile, extensions of such a one-dimensional TASEP are of great 
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interest, since they may be applied to real systems such as protein synthesis [13], bio- 
molecular motors [151 US] traffic flow [T3] , and surface growth [T7] . 

More recently, there are two studies involving a special generalization of the open 
TASEP. Here, the reservoir, or "pool," from which particles are injected into the lattice 
is finite. In particular, the total number of particles on the lattice, iV, plus the number 
in the pool, Np, is a fixed constant: Ntot- Such a constraint can be understood in 
the context of protein synthesis as having a finite number of ribosomes (which model 
particles in TASEP) in a cell[l8], or in the context of traffic (with particles) as 

the "parking garage problem" [T9|. The two studies differ in how particles are moved 
from the reservoir into the lattice, resulting in quite different behavior. In particular, 
the former investigated the properties of N and the average particle current, J, when 
Ntot is varied. While J (Ntot) displays no surprises, N (Ntot) shows some remarkable 
features (when the entry/exit rates are chosen so that, in the Ntot limit, the 

system is in a high-density phase or the "shock" phase). Using simple arguments of 
self-consistency, as well as more sophisticated domain wall (DW) theory [HI [9], many of 
these phenomena can be reproduced [18] . In this paper, we extend this investigation in 
significant ways and provide better insight into the effects of imposing a fixed Ntot 
on TASEP. In particular, we show that, although the previous theory for (Ntot) 
appears quite adequate everywhere, the agreement is deceptive in certain cases. Our 
improvements are not merely incremental; they provide a full understanding of the 
phenomenon of "shock localization." As a consequence, our prediction of the average 
density profile, which is entirely different from that of the simple DW theory, is in 
excellent agreement with simulation results. 

We should note that shock localization has been observed previously [15j. However, 
the underlying mechanisms are distinct. In particular, the earlier studies focus on 
TASEP's with Langmuir kinetics, i.e., one with no particle conservation. To model bio- 
moluecular motors, which can attach and detach from a microtubule, it is natural to let 
particles appear and disappear with various rates along the entire lattice. As a result, 
in approaches that use DW theory and shock dynamics [16], the shock is understood to 
be localized by a subtle interplay between adsorption and desorption of the particles. 
By contrast, our study here involves a TASEP with particle conservation, not only 
within the lattice, but also including the reservoir. The shock is localized through a 
much simpler mechanism: the interplay between the lattice and the reservoir. As will 
be shown below, the mathematical methods used are quite similiar: site dependent 
hopping rates for the domain wall. 

This paper is organized as follows. In the next section, we provide the details of the 
model and summarize previous results. Section [3] is devoted to both new Monte Carlo 
data and the improved DW theory. A summary and outlook for future studies form the 
concluding section. 
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2. Model definition and previous findings 

The standard open TASEP consists of a one-dimensional lattice of sites labeled by 
z = 1, ...,L., each of which can be vacant or occupied by a single particle. Thus, the 
configurations can be specified by the set of occupation numbers {ui}, with = 0, 1 
being the particle content of site i. Thus, N = T^iUi. The (random sequential) dynamics 
is implemented by choosing a particle at random and moving it to the next site with 
unit rate, provided the target site is not occupied. For a particle at site L (the right 
boundary), it leaves the system with rate (3. At the left boundary, a particle can enter 
the system at site 1 with rate a. When the system reaches a stationary state, the 
average overall density 

P = {N) /L 

will be constant. As we vary a and the system can be found in three different phases: 
high density (HD), low density (LD), and maximal current (MC). The MC phase prevails 
if both a and j3 are > 1/2. In the thermodynamic limit, p is 1/2 and J, the average 
current, is 1/4, regardless of (a,/3). If /? < 1/2 and a > /3, the system settles into a 
HD phase, with p = 1 — (3 and J = — (3). Due to particle-hole symmetry, the LD 
phase is similar, with p = a and J = a{l — a), for a < 1/2 and (3 > a. The transitions 
across the phase boundaries HD-MC and LD-MC are continuous. On the line a = (3, a. 
HD region coexists with a LD one, separated by a microscopic interface, known as the 
shock (a term aptly describing a car driving from a high-speed, low-density region into 
a traffic jam). The position of the shock wanders due to fluctuations, so that the long- 
time average of p is again 1/2. Systems displaying such co-existence are often referred 
to as being in the "shock phase" (SP). A concise summary for standard TASEP is 

f 1/2 MC;SP a,(3> 1/2; a = (3 < 1/2 
Poo(a,/3) = < a LD a < min (/3, 1/2) (1) 

[1-/3 HD /? < min (a, 1/2) 

where we have used the subscript oo to remind the reader of its relationship to our 
constrained system {Ntot — ^ oo). The expression for the average current is, in all cases. 

Poo (1 - Poo)- 

In the constrained TASEP, the entry-rate is chosen to depend on Np, the number of 
particles in the reservoir. Denoting this rate by «£//, we characterize such a dependence 
through a function /: 

a,ff = af{Np). (2) 

The advantage of this form is that, by choosing / ^ 1 in the limit of large argument, 
we will recover the standard TASEP with the parameters {a, (3) when the total number 
of particles in the system 

Ntot = Np + N (3) 

is unlimited. In [19j, the reservoir models a "parking garage" and the lattice models 
a road, so that / is chosen to be the simplest function with the desired asymptotic 
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property: / {Np > 0) = 1. Of course, / (0) = 0, since no car can emerge from an empty 
garage. By contrast, the motivation in [18] is the modehng of initiation (a ribosome 
attaching onto an mRNA, to begin the process of protein synthesis) that is hmited by 
the ribosome concentration in the cell. Thus, they chose / (x) oc x, for small x. Much 
of that study was based on the specific function 

/(iVp)=tanh(iVp/iV*) (4) 

where A^* models some cross-over level. As a simple starting point, it was chosen to be 
the average number of particles in the standard TASEP, i.e., A^* = poo L. Thus, 

{L/2 > 1/2; a = /3 < 1/2 

aL a < min(/3,l/2) (5) 

{1-I3)L (3 < min {a, 1/2) 

associated with the MC, LD, and HD states, respectively. As a result, it is more 
convenient to express this function in terms of an intensive control parameter 

Ptot = Ntot/L (6) 

and regard / as a function of p. For example, for the LD case, we have a^/f = 
atanh [(ptot — p) /«]• The main focus in [18j are p {ptot', Ci, P) and J {ptot', en, P), with 
the knowledge that p(oo;a,/9) = poo (a,/5) of Eqn. ([T]). Since predictions for J follow 
those for p, it is sufficient to focus only on the latter. In any case, the variations in J 
are relatively minimal and not as spectacular as those in p. 

To appreciate the different phenomena displayed by the constrained TASEP, it is 
important to keep in mind that the effective entry-rate varies from to a, as ptot is 
increased from to oo. Thus, if a -C min 1/2) (i.e., the "LD case"), then phase 
boundaries are neither traversed nor approached asymptotically. On the other hand, if 
the {a, (3) are chosen so that the unconstrained TASEP is deep in the HD or MC phase, 
we will cross a phase boundary, so that p (ptot) is expected to display two branches. 
In the "HD case", there is a third branch, which reflects the presence of coexistence, 
similar to the dependence of pressure on specific volume, for a liquid-vapour system 
below the critical point. In this branch, cte// = (3 while p is essentially linear in ptot'- 
p = aptot/ (1 — /9 + Oi). In general, the finite size effects, by the time L reaches 1000, are 
hardly noticeable. Using simple self-consistency arguments, the predicted p {ptot) is in 
reasonably good agreement with simulation data [H]. For example, setting p to Oe// for 
the "LD case" (in accordance with Eqn. ([I])), we have p = atanh [{ptot — p) /«]• This is 
a transcendental equation for p, much like the one for magnetisation as a function of the 
external field in the "zeroth" approximation of the Ising model, and leads to p {ptou «)■ 
The most challenging case is "SP," both for simulations and analytic understanding. 
First, the system takes a long time to settle into stationary states. Second, p {ptot) 
shows two cross-overs before saturating at 1/2. Furthermore, the location of the second 
cross-over depends strongly on L. Meanwhile, DW theory for the standard TASEP 
^ accounts for some effects of finite L and provides a more detailed prediction than 
Eqn. ([T]). We denote this result of DW theory by pDw {ci,P,L). Replacing a in this 
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formula by a^jf = a tanh [2 {ptot — p)], another implicit equation for p is established 
[18]. Solving such an equation leads to a p (ptot) that contains all the surprising features 
observed. Indeed, the predicted p (ptot) agrees very well with the data {a = P = 0.25 and 
L = 1000) everywhere, except in the middle of the second cross-over, where the worst 
disagreement is about 6% [H]. Unpublished data for a = /5 = 0.25 and L = 100 and 
analytic results of p {ptot) showed [20j that this second cross-over essentially vanished, 
while the maximum disagreement there between theory and simulations is about 7%. 
In Fig. [H the solid orange symbols represent these two sets of data, while the two lines 
are predictions - with no fitting parameters - from this "simple" domain wall (SDW) 
theory. Given such remarkably good agreements, it may seem an overkill to pursue 
this problem further. However, as we will see in the next section, the agreement in 
[T8] is deceptively good. Had a more sensitive quantity like the average density profile 
(pi = (ui)) been considered, a glaring discrepancy would have been exposed. 
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Figure 1. Average overall density as a function of Ntot for L=100 and L=1000 with 
a = /3 = 0.25. 



3. Domain wall theory for constrained TASEP 

The domain wall (DW) theory used in \181 is formulated for a constant entry (and exit) 
rate [9] , appropriate for the standard, unconstrained TASEP. Referring the reader to the 
details in [9], let us summarize the key points here. The central premise of this theory 
is to let the configurations of the system be approximated by two regions of low/high 
densities, separated by an interface (i.e., DW) of zero "intrinsic width." To be specific, 
a region of low (local) density, p_ (= a), on sites up to k, is connected to a region of 
high density, p+ (= 1 — /5), on sites e [k + l,L], so that only a single integer (fc) is 
used to label each configuration. The dynamics of the DW is implemented through its 
drifting rates to the right/left, D±, dictated by both the densities on either side of the 
shock and the particle fluxes from the two open ends, j±. The final result is a master 
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equation for the probability to find the DW at site k and time t, P {k, t) : 

dtP {k, t) = D+P {k - 1, t) + D_P {k + 1, t) - p+ + D_) P {k, t) (7) 
in the bulk {k G [l,L - 1]), with 

p+ — P- 1 — p — a 

D_ = = ^(^-^) (9) 

— 1 — j3 — a 

When the DW gets to the boundaries, it reflects back into the system, so that 

dtPiO,t) =D^P{l,t)-D+P{0,t) (10) 
dtP {L, t) = D+P {L - 1, t) - D_P (L, t) (11) 

and P {k < 0,t) = P {k > L,t) = 0. The stationary solution is a simple exponential: 

P* (k) = r-^/U (12) 

where 

D_ _ a {I -a) 

and A/" is a normalization factor. It is clear that, for the HD/LD phases, r is 
greater/lesser than unity and the DW is localized at the left/right end of the lattice. 
On the other hand, r is unity for SP, so that the DW can be found anywhere on the 
lattice. 

Turning to our problem of a constrained TASEP, we recognize that, especially if ptot 
is not large, the feedback from the lattice occupation cannot be neglected. As indicated 
above, this is incorporated in the previous study [IB] by introducing an a^ffi^p) and 
solving a self consistent equation for p. But this p is just the average overall density, so 
that the changes in a^jf during a simulation run are completely unaccounted for. Here, 
we wish to investigate how important are the density fluctuations, given that they feed 
back into fte//- To distinguish the fluctuating entry rate from the constant in the SDW 
theory, let us denote the latter by af//'^, so that we have explicitly, 

PsDW (k) oc r;/^ ; r^ff = af/}^ (l - af/J^) / /3 (1 - /3) . (14) 

Now, to account for how this feedback affects the diffusion of the shock in full is 
non-trivial. However, the essentials of the physics are clear: The feedback stabilizes 
the DW, even in the vicinity of SP: If the DW wanders too far to the right (larger 
k), there will be more particles in the pool, resulting in a higher entry-rate. This 
in turn enhances the flux j_ and drives the DW to the left (smaller k). A similar 
stabilizing action occurs if the DW wanders too far to the left. To implement this 
idea and formulate a "generalized" domain wall (GDW) theory, we will make a drastic 
approximation (neglecting time delays, etc.), consisting of three ingredients: 

(i) replacing a by Og// everywhere 

(ii) substituting p = p- (k/L) + p-|_ (1 — k/L) into Oe// (p) 
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(iii) letting p_ be Oe// and obtaining a fc-dependent ae//,fc through a self-consistent 
equation. 

The last point is subtle and deserves clarification. As a^jf = af [Np), i.e., 
atanh [{ptot — p) /«] here, we have an imphcit equation for determining ae//,fc : 

aeff = ataxi\i[{ptot- P-{k/L) - p+{l-k/L)) /a] (15) 
aeff,k = atanh [{ptot - a^jf^k (k/L) - (1 - /5) (1 - k/L)) /a] (16) 
Once aeffk is found, it will enter in defining the fc-dependent drift rates: 

D.^^^J^^ (17) 

1 - /? - aeff,k 

The result of these considerations is a generalized master equation for P {k, t) which, 
in the "bulk," reads 

dtP{k) = D+,k-iP{k - 1) + D-,fe+iP(A; + 1) - (/}+,, + D^,k)P{k) (19) 

(with the variable t suppressed). For the boundary equations, there is a further 
complication which we must account for. When the resources are so limited that the 
lattice cannot be fully filled at density p+, i.e., when Ntot < P+L, the DW cannot be 
located arbitrarily far to the left. Thus, the interval available to k is [/cmin,-^], where 

Ntot 

kmin = L ~ —. (20) 

Note that k = fcmin corresponds to the configuration with A'p = 0, and so, p_ = 
ae//(0) = 0. With this limitation in mind, we obtain the boundary terms for the 
master equation: 

dtP{krnin) = D-,fc^,„+lP(fc^in + 1) " Z^+,fc^,„P(A;^in) (21) 

dtP{L) = D+,L-iP{L - 1) - D^,lP{L) (22) 

Despite the extra complications, the stationary distribution can be found 
analytically |16j. Since the configuration space is one-dimensional, we have, in general, 
a recursion relation P*{k) = [D^^k+i/D+,k] P*{k + 1). Thus, we arrive at the solution 
for our constrained TASEP: 



^*(^)-n^ (23) 

l=k 

and, with proper normalization. 



k-l 



\m=k,^in e=m+l ^' m=k+ie=k+l ' / 

Since this P*{k) carries information on whether the DW is located before or after a 
given site i, we can compute the average density profile: 

i L 

p,= J] (1 - P)P*{k) + a,ff,kP*{k) (25) 
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and of course, the overall density is given by p = ^pj/L. Needless to say, these results 
are quite different from those in the simple DW theory, e.g., Eqn. f|T2|) . 

To highlight how differently the two predictions compare with simulation data, we 
choose just one point in parameter space - a point in the "third branch" of p {ptot) in 
the HD case. Here, the essential physics revealed by simulations is that the number 
of particles in the pool remains more or less constant as Ntot is increased, with the 
extra particles being absorbed by the lattice to form a high density region near the 
exit. Furthermore, the shock on the lattice is seen to be quite localized, as revealed 
by the average density profile. Fig. [2] shows this behavior for the case of L = 1000, 
Ntot = 800, a = 0.75, and (3 = 0.25. In this figure, we also plot the predictions from 
the two theories, Eqns. (11211241) . Though both provide quite good agreement with 
the overall density (area under the profile), it is clear that shock localization cannot 
be achieved by using an exponential P*{k). By contrast, accounting for the feedback 
through site-dependence drift rates for the DW|16]. the shock is localized. Given how 
crude our approximations are, it is remarkable that the agreement with data is so good. 
We conclude that even this simple-minded level of accounting for feedback can capture 
the essence of shock localization. 




Figure 2. Average density profile for L — 1000 witli Ntot — 800, a — 0.75, and 
/? = 0.25. 

Apart from the dramatic improvements to the profile here, our GDW theory also 
provided better fits to both p (ptot) and the profile in the SP case. Since our theory 
accounts for some feedback, we are not surprised that it is more successful at dealing 
with the large fiuctuations associated with the SP case. Carrying out the computations 
detailed above, we find results that are in surprisingly good agreement with simulation 
data, as the open circles in Figs. [1] and [3] show. While we show two cases in the former 
(L = 1000 and L = 100), the latter contains only the profile for one case - ptot = 2 
- in the L = 1000 samples. Although the improvements over SDW here are not as 
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Figure 3. Average density profile for L = 1000 witli Ntot = 2000, a = /3 = 0.25. 

Another useful perspective on the difference between the SDW and GDW theories 
is the following. Consider the gradient of the density profile (which is intimately related 
to P*{k), of course). For a discrete lattice, this is 

Apfe = Pk- Pk-i ■ 

In SDW, it is simply (1 — /3 — afjj'^)P|^(4/(fc), and so, is a pure exponential, since c^eff^ 
is a constant. In other words. In Ap^ is linear in k (with coefficient Inrg//). By contrast, 
in GDW, an extra k dependence appears through cte/f^k in (1 — /5 — aeff^k)P*{k). From 
Fig. [2], it is clear that the non-linear terms are significant and Ap^ is close to a Gaussian. 
On the other hand, it is also clear that Fig. [3] shows that InAp^ is essentially linear. 
That GDW can account for such sharp differences is remarkable. At a deeper level, 
to understand how a simple k dependence in ae//,fc can be so successful will require a 
thorough analysis of the detailed properties of / (iVp). Beyond the scope of this paper, 
this study is being undertaken and will be reported elsewhere. 

4. Summary and outlook 

In this paper, we revisit the constrained totally asymmetric simple exclusion process 
proposed in [18j. Instead of particles entering/exiting the lattice from/to an infinite 
pool, this system has a fixed and finite total number of particles, Ntot ^ the sum of 
those on the lattice (A^) and in the pool (Ap). Further, the entry rate cte// is assumed 
to be a simple function of Ap, in such a way that a^jf oc Np if the pool has few 
particles, yet saturating at a constant a for Ap ^ oo. The main interest is how A^ 
varies as Ntot is increased, when the exit rate is always fixed at j3. In [18], Monte Carlo 
simulation data showed that N {Ntot', C(, j3) displays a variety of properties, depending 
on the parameters {a, (3). A simple domain wall theory was shown to provide good 
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predictions {zero adjustable parameters) for all the data. However, there were lingering 
doubts for the rationale behind such a simple theory. Here, we undertake to improve 
upon these theoretical considerations. 

By taking into account the fluctuating feedback into the entry rates from the total 
particle occupation on the lattice, we formulated a "generalized" domain wall theory, 
in which the hopping rates for the domain wall depend on the location of the wall 
(similar to the spirit, but differing in the physics of earlier studies |16]). Although 
this approach is still somewhat heuristic, its predictions for N {Ntot', en, (3) are better 
than those in [18]. Furthermore, we find that the agreement in [18] turns out to be 
deceptively good for certain cases. In such situations, the previous predictions for a more 
detailed quantity - the density profile - can be seriously erroneous. By contrast, the 
generalized domain wall theory provides excellent results for the profiles in all cases. In 
physically understandable terms, the feedback mechanism serves to localize the domain 
wall. While this feedback plays a small role in most cases, its effects are all-important 
in other cases, with the most dramatic example displayed in Fig. O Our conclusion 
is that even a naive inclusion of the feedback improves significantly our understanding 
of the fluctuations in a constrained TASEP. Of course, further progress can be made, 
such as a more systematic approach to accounting for all the fluctuations in this system. 
In particular, the quantitative aspects of this feedback are clearly associated with the 
details of /, specifically, its derivatives. Beyond static properties, there are undoubtedly 
many interesting dynamic phenomena yet to be discovered. For example, we are aware 
of non-trivial effects [21] induced by the constraint on the power spectrum associated 
with the time trace of [10] - a phenomenon that may be understood through an 
appropriate extension of this theory. 

Finally, we should remark that one of the motivations for studying TASEP with 
finite resources comes from the potential applications to protein synthesis in a cell. In 
a real biological system, the supply of ribosomes (which are modeled by particles in 
TASEP) is finite. In case this supply is low, there may be observable consequences for 
translation. For that application, the above TASEP needs to be generalized, to include 
exclusion at a distance (particles covering more than one site) and inhomogeneous hop- 
ping rates (associated with a non-trivial sequence of codons) [13]. Furthermore, there 
are many different genes, as well as many copies of each. Thus, we would face a system 
of multiple copies of TASEP's, with different lengths and hopping rates. How these 
systems are affected by finite resources will pose many interesting new challenges, es- 
pecially on the theory front. It is clear that our study here is but a very small step 
towards the understanding of protein synthesis in vivo. 
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